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ABSTRACT 

We study the aggregation of peptides using the discrete molecular dynamics simulations. 
At temperatures above the a-helix melting temperature of a single peptide, the model pep- 
tides aggregate into a multi-layer parallel /5-sheet structure. This structure has an inter- 
strand distance of 4.8A and an inter-sheet distance of lOA, which agree with experimental 
observations. In this model, the hydrogen bond interactions give rise to the inter-strand 
spacing in /3-sheets, while the Go interactions among side chains make /3-strands parallel to 
each other and allow /5-sheets to pack into layers. The aggregates also contain free edges 
which may allow for further aggregation of model peptides to form elongated fibrils. 



I. INTRODUCTION 

Protein misfolding and polypeptide aggregation are in the focus of interdisciplinary sta- 
tistical physics because of their relevance to amyloid diseases such as Alzheimer's disease, 
Parkinson's disease and Huntington's disease. Even though polypeptides related to these 
diseases share no sequence or secondary structure similarity, they can aggregate into insolu- 



ipoA in diameter, 
reveal the com- 



ble fibrils which share some structural features. These fibrils are typical 

n 

and several thousand Angstroms in length [JJ]. X-ray diffraction studies 
mon structural features for these amyloid fibrils: the presence of a 4.7-4.8 A inter-strand 
spacing along the fibril axis and a 9-10 A inter-sheet spacing perpendicular to the fibril 
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axis |^|5|. Although advances have been made toward understanding the structural charac- 
teristics of the fibrils and the mechanism of fibril formation, our knowledge of the detailed 
fibrillar structure and mechanisms of amyloid assembly is limited. 

Molecular dynamics provides a way to study the aggregation mechanism at the molecu- 
lar level. However, continuous all-atom molecular dynamics simulations with realistic force 
fields in a physiological solution are not fast enough to monitor a complete aggregation pro- 
cess from monomers to fully-formed fibrils. Recently, a discrete molecular dynamics (DMD) 
algorithm 0, Q| using a coarse-grained peptide model has been successfully implemented 
to study protein folding thermodynamics and protein folding kinetics j^. This computa- 
tionally fast and dynamically realistic simulation technique has also been applied to study 
the aggregation of a small number of Src SH3 domain proteins P] and the competition of 
refolding and aggregation of four- helix bundles jlo| . 

Here we study the aggre gatio n of a large number of peptides. We choose 40-amino acid 
amyloid [3 peptide (A/?i_^n protein data bank (PDB) [3] access code 1BA4), which 
is associated with Alzheimer's disease, to construct model peptides. Our results show that 
model peptides can aggregate into multi-layer /9-sheet structures with free edges which 
may enable further fibrillar elongation. The computed diffraction pattern of our simulated 



multi-layer /3-sheet is consistent with experimental observations 
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3. 



II. METHODS 



A. Two-bead model 



1. Geometry of model peptide: beads and permanent bonds 

We model each amino acid in the A/5i_4o peptide by two beads - bead representing 
backbone atoms and Cp bead representing side chain atoms (for glycine, Cp is absent). Each 
bead has an index i indicating the position of amino acid in the sequence starting from the 
N-terminus. The geometry of the peptide is modeled by applying permanent bonds among 
these beads These bonds include covalent bonds between C^i and Cpi, peptide bonds 
between Cai and Ca{i±i), additional constraints between Cpi and CQ,(j-|-i), and also between 
Cai and Ca(i±2) (Fig. H}. These additional constraints are introduced to model angular 
constraints between side chains and the backbone. 
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All permanent bonds are realized by infinitely high potential well interactions between 
the related beads 



Dij{l - aij) < \ri - r,| < 0^(1 + 
-oo otherwise. 



ybond ^ ^ " ^ijj ^ \'i 'J\ I ^ijj , 



Here Dij is the bond length between beads i and j, and is the relative deviation of this 
bond length. The average lengths for these bonds can be obtained from statistical analysis 



of distances within the y4/3i_4n_NMR structures 



and their relative deviations llTll used in our model. 



Table |l] presents the average lengths 



2. Interactions between beads: Go model 

Typically the Go potentials are used to model proteins with well-defined globular 

native states. Side chains which form contacts in the native state (native contacts) experience 
attractive Go potential. However, y4/5i_4o peptide is "natively unfolded". NMR studies 
suggest that in hydrophobic environments the APi_4q peptide assumes mostly a-helical 
conformation Fig. 01 (a) shows one of these NMR structures. Therefore, we apply Go 
potentials to preserve this well-defined, mostly a-helical structure of the v4/3i_4o peptide. 
In our two-bead model a native contact is defined when two Cjj beads are closer than 
jjGo _ 7 within the NMR structure of the A(3i-4o peptide. All the C/3 beads can not 
be closer than the hard-core distance -Dj^^ = 4.5A. In particular, the structure-specific Go 
potentials make the side chains indexed by i within the a-helix region of v4/?i_4o peptide 
attract side chains z ± 2, z ± 3 and z ± 4. Fig. El (b) shows the native contact map for the 
NMR structure of /l/3i_4o peptide shown in Fig. El (a). 

To study the aggregation we need to simulate also the interactions between different 
peptides. We apply Go potentials for Cp beads in different peptides by an assumption that 
two amino acids which interact with each other in a single peptide will interact in the same 
way in different peptides. For example, amino acids 16 and 19 form a native contact in 
the NMR structure. Thus, amino acids at 16 and 19 of peptide 1 will experience attractive 
Go-type interaction with amino acids 19 and 16 of the peptide 2, respectively. The strength 
of Go interactions is set to unity e*^" = 1. 



3 



3. Interactions between Ca heads: hydrogen bond 

For many globular proteins it has been observed that the number of backbone hydrogen 
bonds for each amino acid does not exceed two Also, whenever two hydrogen bonds are 
formed in a particular peptide block they are approximately parallel to each other. In order 
to incorporate these two facts in our model we introduce two criteria for hydrogen bond 
formation: (i) that each Ca bead can form up to two effective hydrogen bonds, and (ii) that 
the two hydrogen bonds formed by the same Ca bead must be approximately parallel. 

We set the hydrogen bond interaction range between two beads to D^^ = 5.0A, and 
their hard-core distance to -D|f^ = 4.0A. We use the following procedure in order to satisfy 
the criteria for the hydrogen bond formation: when two beads, A and B, come to a 
distance D^^, we check for any existing hydrogen-bond partners of A and B. If both beads 
A and B have no existing hydrogen partners they can form a hydrogen bond automatically. 
If one of the beads, for example A, already has one partner, Ai, and the distance between 
the bead Ai and the bead B is within the range of 8.7-lOA (i.e. the angle between vectors 
AAi and A]3 is within the range of 120°-180°), the bead A can form another hydrogen 
bond with bead B provided that either the bead B has no existing hydrogen bonds or its 
single hydrogen bond partner, Bi, has a distance with bead A in the range of 8.7-loA (see 
Fig. 121). If one of beads A and B or both already have two hydrogen bond partners, the 
pair will proceed with a hard-core collision without forming a new hydrogen bond. When 
a new hydrogen bond is formed between beads A and B, new hydrogen bond partners are 
recorded for these two beads, and whenever a bead gets two hydrogen bond partners an 
auxiliary bond is formed between these two partners. Every auxiliary bond can fluctuate 
within the range of 8.7-lOA to keep two hydrogen bonds within the angle 120°-180° and 
it cannot be broken unless one of the two hydrogen bonds is broken. A hydrogen bond 
between beads A and B can be broken when these two beads move away from each other to 
a distance of D^^ and their kinetic energies are higher than e^^ . When a hydrogen bond 
is formed or broken, the velocities of the beads A and B change in order to conserve energy 
and momentum, such that their kinetic energy increases or decreases by the value e^^ . We 
set e^^ = 3 as it was chosen in Ref 19| for Src SH3 domain. 
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B. Computed diffraction 



For the typical conformation of A/3i_4o peptide aggregation structure, we calculate the 
intensity of diffraction pattern using the elastic diffraction formula P| in order to compare 
with experimental results ll^ • 

I{kf) = \Y^exp{z{kf-k,)-r,)\^ (2) 
j 

where ki is the wave vector of the incidental X-ray, kf is the wave vector of the outgoing 
X-ray, fj is the position vector of jth bead, and the summation is taken over all the and 
Cj3 beads in the structure. 

We chose x axis perpendicular to the /3-sheets, and y axis along the fibrillar axis which 
is perpendicular to the /^-strands in the /5-sheets(Fig. [7| (a)). The incoming X-ray with lA 
wavelength goes along z axis and the diffraction pattern is collected on a. x — y plane behind 
the aggregate sample. The deflecting angle, 6 = cos~^{kf ■ ki/k"^), ranges from 0.05 to 0.25 
in radians in order to detect the periodicity of 4A to 20A in the aggregate structure. Since 
amyloid fibrils consist of bundles of /3-sheet chains which are twisted along the y-axis, there 
is no preferred orientation in the x — z plane in the X-ray diffraction experiments. We rotate 
the structure candidate around the y (fibrillar) axis n times by angle 2-7r/n and add all the 
diffraction intensities to obtain a final pattern. We take n = 20 in the present study. 



III. RESULTS FOR A SINGLE PEPTIDE 



As an initial test of our model peptide, we perform DMD simulations of a single pep- 
tide to test whether a peptide with random coil conformation recovers the observed NMR 
structure. The model peptide is slowly cooled from Tj=1.00, which is high enough to ren- 
der the peptide as a random coil, to different target temperatures Tj=0.60, 0.55, 0.25. 
For each target temperature we make 10 trials starting with different initial conformations. 
When Tt<0.40, the segment Q15-V36 adopts an a-helix or two pieces of left-handed and 
right-handed a-helices. This artifact is observed because our simplified two-bead model does 
not distinguish between different handedness. At Tt=0.40, the N terminus adopts mostly a 
random coil conformation. As Tt<0.35, the model peptide starts to approach to its ground 
state which is an a-helix with a single handedness along the entire peptide chain. Therefore, 
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as expected within a certain temperature range around T=0.40 during the coohng process 
the model peptide adopts partial a-helical conformation similar to the observed one in NMR 
experiments. 

We also study the equilibrium behavior of a single model peptide at different temperatures 
by measuring the specific heat as a function of temperature. At each sampled temperature 
we start with a ground state conformation and perform DMD for 10^ simulation time steps 
to equilibrate the system, followed by additional 10'' time steps for the calculations. Figs. E] 
(a) and (b) show the potential energy and specific heat as a function of temperature for a 
single model peptide, respectively. The melting of a-helix is non-cooperative which can be 
concluded from the broad peak between and T^^O.SS in the specific heat curve 

(Fig. 111(b)). T/v corresponds to the structural transition from an a-helix to a random coil 
for the first 14 amino acids starting from the N-terminus, while Tm corresponds to the 
melting of the a-helix in the segment Q15-V36. is higher than T/v because there are 
more attractions among the side chains in the segment Q15-V36. 

IV. RESULTS FOR MULTIPLE PEPTIDES 

In the study of aggregation of many identical peptides, we perform simulations of 28 
peptides in a cubic box with the edge of 200A and periodic boundary conditions. Initially, 
all the peptides are placed on a grid and randomly oriented (see Fig. El (a)). Then we 
equilibrate the system at various temperatures: Tj-=0.4, 0.5, 1.20. 

At temperatures lower than the melting temperature of a single peptide, peptides 
in our model aggregate into amorphous structures where individual peptides preserve most 
of their a-helical segments as in Fig. El (b). When the temperature is higher than T^,, 
peptides start to aggregate into more ordered structures. When the temperature is higher 
than 1.10, there is no stable aggregate (this threshold temperature depends on the peptide 
concentration). At a temperature range between 0.55 and 1.10, the model peptides can 
aggregate into multi-layer /3-sheet structures. Figs. IHl (a) and (b) show the time evolution 
of the conformation obtained from DMD simulation at temperature 0.90. In Fig. [71 (a) we 
illustrate the setup of diffraction computation and in Fig. |3 (b) we present the calculated 
diffraction pattern. The relative sharp and intense 4.8 A meridional reffection corresponds 
to the periodic packing of /3-strands along the fibril axis, and the weaker lOA equatorial 
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reflection corresponds to the distance between /3-sheets. In Fig. |7| (c) we show the calculated 
pair correlation function for the same /3-sheet structure. The peaks around 4.8A and lOA 
correspond to the average inter-strand and inter-sheet spacings, respectively. 

In order to study the thermostability of this 3-layer /3-sheet structure, we slowly increase 
the temperature to T=2.0 which is higher enough to melt the aggregate. Figs. |Hl(a) and (b) 
show the time evolution of temperature and the temperature dependence of potential energy 
of the system during melting and dissociation of the /3-sheet structure, respectively. As 
temperature increases from 0.90, the aggregate becomes less stable. At temperature around 
T=1.15 ± 0.05, aggregate starts to dissociate. At temperatures higher than Td=1.20 ± 0.05 
the dissociation completes. 

If we assume that the temperature 0.9 at which the aggregation of /3-sheet is observed 
corresponds to physiological temperature, 310 K. At this temperature our single model 
peptide exists in a random coil conformation, which corresponds to experimental observa- 
tions [20, |2l| that in aqueous solution at physiological temperatures /l/3i_4o peptides adopt 
mostly /3-sheet and coil conformations. The temperature of the /3-sheet dissociation, 1.2, 
corresponds to 413 K. Temperature T=0.40 at which our model peptide acquires a-helical 
conformation corresponds to very low physical temperatures which can not be observed 
experimentally. 

V. DISCUSSION AND CONCLUSION 

In the test of our coarse-grained model of v4/3i_4o peptide, we find that the model peptide 
most resembles the NMR structure of yl/3i_4o peptide around T=0.40. The existence of an 



optimal temperature range 
other coarse-grained models 



23 and 



br protein refolding is also observed in experiments 
2^. Below T=0.40 the N-terminal region of our model peptide 
mostly adopts an a-helical conformation. However, in the present study of aggregation we 
are focused on temperatures above 0.40 as the peptides are generally partially or completely 
unfolded to initiate the aggregation process j2J]. 

In studies of multiple peptides, we demonstrate that peptides aggregate into amorphous 
structures (Fig.El(b)) around T=0. 50 or multiple-layer /3-sheet structures (Fig.|ni(b)) around 
T=0.90. In the amorphous structures, individual peptides tend to preserve most of the a- 
helical structure, while in the /3-sheet structures the /3-strands tend to be parallel. Since the 
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Go interaction for an a-helix favors the formation of contacts between amino acids i and 
z ± 3, the aggregates with a parallel alignment have lower potential energies. 

There are caveats due to the simplicity of the two-bead model used in our study. Each 
amino acid is represented by only two beads, which do not allow for an accurate description 
of the backbone. The backbone in this model is too flexible, which introduces some artifacts 
into conformations of aggregates composed of small number of peptides at low temperatures, 
such as dimers, trimers and tetramers. 

An additional problem is that the chiral symmetry of each amino acid is not considered in 
this model. As a result, we observe two a-helices with opposite handedness. As the /l/3i_4o 
NMR structure consists of two a-helices and a hinge in between, there are four low energy 
states with combinations of different handedness within the region of a-helices at T=0.40. 
The conformations with mixed handedness appear with lower probabilities since they have 
higher potential energies due to the loss of native contacts and hydrogen bonds in between 
the two a-helices of different handedness. 

Also, due to the simplicity of the two-bead model, we do not account for specific struc- 
tural features of APi^^o peptides, such as the salt bridge between D23 and K28 For the 
same reason, we can not expect to explain the differences in aggregation pathways between 



/l/3i_4o and /l/3i_42alloforms [26|, nor study subtle aggregation differences due to amino 
acid substitutions [27 1. We show that the DMD algorithm using a simplified peptide model 
can reproduce the formation of /5-sheet structures of 28 peptides with free edges for further 
fibrillization. Our study shows that it is possible to investigate in detail the aggregation of 
several dozens of peptides using DMD simulations and the coarse-grained model for peptide 
structure. Both the number of peptides and the complexity of the model j^l can be signifi- 
cantly increased within realistic computational constraints. Thus we regard this study as a 
first step toward developing a realistic model of Af3 peptide aggregation. 
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TABLE I: Permanet bonds in the two-bead model 
Bond Bond length (A) Deviation (%) 

Co^i - Cp, 1.55 2.4 

3.82 3.1 
C/3i - C'«(i±i) 4.66 6.5 

Cai - C„(i±2) 5.65 14.8 




FIG. 1: Schematic diagram of two-bead model. Each amino acid in the ^/9i_4o peptide is rep- 
resented by two beads: Ca bead represents backbone atoms and Cp bead represents side chain 
atoms {Cp is absent for amino acid glycine). The geometry of the peptide is modeled by applying 
permanent bonds among these beads: covalent bonds (bold lines), peptide bonds (thin lines) and 
additional constraints (dashed and dotted lines). Interactions between side chains are modeled 
by Go potentials between Cp beads, and interactions between backbone atoms are modeled by 
hydrogen bond interactions between Ca beads. 
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A 



B 



FIG. 2: Model of a hydrogen bond. Existing hydrogen bonds AAi and BBi are shown in bold 
lines. When the beads A and B come to a distance SA, a new hydrogen bond (dotted line) may 
form if the distances AiB and BiA satisfy inequalities 8.7A < AiB < 10. oA and 8.7A < BiA< 
10. oA. If the bond AB is formed, the auxiliary bonds AiB and BiA (dashed lines) are formed 
simultaneously. These bonds can fluctuate within the interval 8.7-loA and cannot be broken unless 
beads A and B move away from each other to a distance SA. If the beads A and B have enough 
kinetic energy to leave the hydrogen bond attraction well, their velocities are changed in order 
to conserve energy and momentum, and the hydrogen bond AB is destroyed simultaneously with 
the auxiliary bonds AiB and BiA. The velocities of Ai and Bi do not change at the moment 
of forming or destroying of hydrogen bond AB. Analogously, if one of the hydrogen bonds, AiA 
or BiB, breaks before hydrogen bond AB, the corresponding auxiliary bonds AiB or BiA also 
breaks. 
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FIG. 3: (a) The NMR structure of A/3i_4o peptide [l]| used to construct the two-bead model 
peptide with Go potentials and hydrogen bond interactions. The picture is created with the 
program Molscript [2^. (b) The contact map for structure (a). Note that the a-helical region is 
from amino acid 15 to 36 (Q15-V36). 
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FIG. 4: Temperature dependence of (a) potential energy and (b) specific heat for a single two-bead 
A/3i_4o rnodel peptide with Go potentials and hydrogen bond interactions. The calculations are 
based on the DMD simulations of 10^ time steps for each sampled temperature. 



(a) 




FIG. 5: DMD simulation of 28 peptides at temperature 0.5: (a) initially, all peptides in an original 
randomly oriented NMR conformation are placed on a grid, (b) An amorphous aggregate obtained 
by DMD simulation at this temperature of 0.5. The simulation shows that most of the a-helical 
segments are preserved during the aggregation at this temperatures. 
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FIG. 6: DMD simulation of 28 peptides at temperature 0.90. Initial conformation is the same as 
Fig. [5] (a) . After five hundred time steps all peptides acquire random coil conformation character- 
istics for T=0.90 (data not shown), (a) Intermediate conformation at temperature 0.90 after 10"* 
DMD simulation time steps, (b) 3-layer parallel /3-sheet structure formed after 2.5 x 10^ DMD 
simulation time steps. This /3-sheet structure contains free edges which may allow for further 
aggregation of model peptides along the y-axis, which is perpendicular to the plane of the figure. 
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FIG. 7: (a) The setup of diffraction pattern computation for the 3-layer /?-sheet aggregate formed 
by 28 peptides shown in Fig. [BJa at a different perspective. The scattering occurs along z-axis 
which is perpendicular to the plane of the figure, (b) Computed diffraction pattern collected on a 
x-y plane behind the aggregate. The pattern is averaged over 20 patterns obtained by successively 
rotation of the aggregate around y-axis by 18°. (c) Pair correlation function for the same aggregate, 
where peaks around 4.8 A and 10 A correspond to inter-strand spacing and inter-sheet spacing, 
respectively. And the peak around 5.7 A is mainly from the correlation between neighboring Cp 
beads. 
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FIG. 8: The melting of the 3-layer /3-sheet structure of 28 Af3 peptides (Fig. El (c)). (a) Time 
evolution of the temperature when the 3-layer /?-sheet is warmed up slowly from temperature 
0.90 to 2.00. The dissociation temperature is Trf=1.20 it 0.05. The insets are the conformations at 
different times/temperatures. Note that the third one shows a completely dissociated conformation 
in 3-d. (b) The temperature dependence of potential energy. 
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